************************************************************
* Project: Defending the Watchdogs                         *
* Authors: Jonathan A. Solis  & Kelebogile Zvobgo          *
* Task: Analyze/Graph Interations                          *
* Date: 7/10/2022                                          *
* Software: Stata MP 17                                    *
************************************************************

*-------------------------------------------------------------------------------
*NOTE: Please contact Jonathan (jsolis@aiddata.wm.edu) if you have any questions
*-------------------------------------------------------------------------------

*******************************************************************************************
*Instruction:
* -This code graphs the interactions found in the manuscript and appendix. 
* -To make it work, simply find the model you wish to run (marked below) and
*  ensure (i) the model and (ii) 'gen byte used=e(sample)' are un-greened/activated
* -To emphasize, only (i) the model and (ii) 'gen byte used=e(sample)' are un-greened/activated
*  and no other models. 
* -Once it is set up, run entire script, and this will generate the figure you want.
* -To get another figure, green the (i) the model and (ii) 'gen byte used=e(sample)'
*  you just used and un-green/activate (i) the model and (ii) 'gen byte used=e(sample)'
*  in which you're interested.
*******************************************************************************************


*##* Load Data
*-------------
clear 
cd "[insert path]"
use DefendingWatchdogs_Dataset
set more off

*##* Set up dataset for Analysis
*-------------------------------

*Special ado file
***search gr0002_3
*-> for lean1 scheme for figures

*Natural log gdppc variables
gen gdppc_ln=ln(gdppc)

*Give Malta a region/colonial backgroun
replace region=5 if country=="Malta"
replace nawest=1 if country=="Malta"
replace colonized=0 if country=="Malta"
replace british=0 if country=="Malta"
replace french=0 if country=="Malta"
replace spanish=0 if country=="Malta"
replace port=0 if country=="Malta"

*Create Interactions*

*Judicial Independence (L&S)
gen lji_poly=lji*polyarchy
gen lji_gold=lji*gold

*Judicial Independence (VDEM)
gen vdem_poly_hc=vdem_ji*polyarchy
gen vdem_poly_const=jconstraint*polyarchy
gen vdem_gold=vdem_ji*gold

*Create Macros*

*DVs (w/ lag)
global mdcensor mdcensor l.mdcensor
global mdharass mdharass l.mdharass
global inetcen   inetcen  l.inetcen

*IV variables (LJI): simple model
global xlist_simple           l.lji l.polyarchy
global xlist_simple_gold      l.lji l.gold 
global xlist_simple_int       l.lji l.polyarchy l.lji_poly
global xlist_simple_gold_int  l.lji l.gold5 l.lji_gold

*IV variables (VDEM-judicial constraint): simple model
global xlist_simple_vdemjc      l.vdem_jconstraint l.polyarchy l.vdemjc_poly
global xlist_simple_gold_vdemjc l.vdem_jconstraint l.gold l.vdemjc_gold

*IV variables (VDEM-hc): simple model
global xlist_simple_vd       l.vdem_ji l.polyarchy l.vdem_poly_hc
global xlist_simple_gold_vd  l.vdem_ji l.gold l.vdem_gold

*IV variables (VDEM-hc-ord): simple model
gen vdem_poly_hco=vdem_jiord*polyarchy
global xlist_simple_vdo       l.vdem_jiord l.polyarchy l.vdem_poly_hco
global xlist_simple_gold_vdo  l.vdem_jiord l.gold l.vdemo_gold

*IV variables: controls
global xlist_ctrls l.gmf_binary l.education l.KOFInGIdj l.gdppc_ln l.aid l.rents_wb l.conflict2 l.hrp_point

*set as panel
xtset ccode year

*set obs to 100000
set obs 100000

*Models 

* * * * * * * * *
*****************
*TRAD CENSORSHIP*
*****************
* * * * * * * * *

*------*
*GLOBAL*
*------*

xtreg mdcensor $xlist_simple_int $xlist_ctrls i.year, fe vce(robust)
gen byte used=e(sample)

*.87

* sum polyarchy if polyarchy>.87&used==1
*
*    Variable |        Obs        Mean    Std. Dev.       Min        Max
*-------------+---------------------------------------------------------
*   polyarchy |        667    .8890135     .013867        .87       .924

* display 667/6410
* .10405616

*-------------------------------*
*East Europe/Former Soviet Union*
*-------------------------------*

*xtreg mdcensor $xlist_simple_int $xlist_ctrls i.year if eeuro==1, fe vce(robust)
*gen byte used=e(sample)

*-----------------------*
*Latin America/Carribean*
*-----------------------*

*xtreg mdcensor $xlist_simple_int $xlist_ctrls i.year if lac==1, fe vce(robust)
*gen byte used=e(sample)

*------------------------*
*Middle East/North Africa*
*------------------------*

*xtreg mdcensor $xlist_simple_int $xlist_ctrls i.year if mena==1, fe vce(robust)
*gen byte used=e(sample)

*------------------*
*Sub-Saharan Africa*
*------------------*

*xtreg mdcensor $xlist_simple_int $xlist_ctrls  i.year if ssa==1, fe vce(robust)
*gen byte used=e(sample)

*--------------------*
*N. America/W. Europe*
*--------------------*

*xtreg mdcensor $xlist_simple_int $xlist_ctrls  i.year if nawest==1, fe vce(robust)
*gen byte used=e(sample)

*----*
*Asia*
*----*

*xtreg mdcensor $xlist_simple_int $xlist_ctrls i.year if asia==1, fe vce(robust)
*gen byte used=e(sample)

* * * * * * * * * * 
********************
**MEDIA HARASSMENT**
********************
* * * * * * * * * * 

*------*
*GLOBAL*
*------*

*xtreg mdharass $xlist_simple_int $xlist_ctrls i.year, fe vce(robust)
*gen byte used=e(sample)

*.54

* sum polyarchy if polyarchy>.54&used==1
*
*    Variable |        Obs        Mean    Std. Dev.       Min        Max
*-------------+---------------------------------------------------------
*   polyarchy |      2,761    .7705444    .1119543        .54       .924

* display 2761/6410
* .43073323

*-------------------------------*
*East Europe/Former Soviet Union*
*-------------------------------*

*xtreg mdharass $xlist_simple_int $xlist_ctrls  i.year if eeuro==1, fe vce(robust)
*gen byte used=e(sample)

*-----------------------*
*Latin America/Carribean*
*-----------------------*

*xtreg mdharass $xlist_simple_int $xlist_ctrls i.year if lac==1, fe vce(robust)
*gen byte used=e(sample)

*------------------------*
*Middle East/North Africa*
*------------------------*

*xtreg mdharass $xlist_simple_int $xlist_ctrls i.year if mena==1, fe vce(robust)
*gen byte used=e(sample)

*------------------*
*Sub-Saharan Africa*
*------------------*

*xtreg mdharass $xlist_simple_int $xlist_ctrls i.year if ssa==1, fe vce(robust)
*gen byte used=e(sample)

*--------------------*
*N. America/W. Europe*
*--------------------*

*xtreg mdharass $xlist_simple_int $xlist_ctrls i.year if nawest==1, fe vce(robust)
*gen byte used=e(sample)

*----*
*Asia*
*----*

*xtreg mdharass $xlist_simple_int $xlist_ctrls i.year, fe vce(robust)
*gen byte used=e(sample)

* * * * * * * * * * * * 
***********************
**INTERNET CENSORSHIP**
***********************
* * * * * * * * * * * * 

*xtreg inetcen $xlist_simple_int $xlist_ctrls i.year, fe vce(robust)
*gen byte used=e(sample)

* * * * * * * 
*************
***V-DEM JI**
*************
* * * * * * * 

*----------*
*Censorship*
*----------*

*xtreg mdcensor $xlist_simple_vdo $xlist_ctrls i.year, fe vce(robust)
*gen byte used=e(sample)

*----------*
*Harassment*
*----------*

*xtreg mdharass $xlist_simple_vd $xlist_ctrls i.year, fe vce(robust)
*gen byte used=e(sample)

*--------*
*Internet*
*--------*

*xtreg inetcen $xlist_simple_vd  $xlist_ctrls i.year, fe vce(robust)
*gen byte used=e(sample)

* * * * 
*******

*------------------------------*
*CREATE MARGINAL EFFECTS FIGURE*
*------------------------------*

matrix b=e(b)
matrix V=e(V)

*see coefficient vector
matrix list b
*see covariance matrix
matrix list V

*save scalar of beta coefficeint X "lji"
scalar b1=b[1,1]
*save scalar of interaction
scalar b3=b[1,3]

*save scalar of variance of X "lji" w/ itself
scalar varb1=V[1,1]
*save scalar of variance of interaction and itslef
scalar varb3=V[3,3]

*make scalar of covariance of X 'lji" and interaction
scalar covb1b3=V[1,3]

*list of what you just created
scalar list b1 b3 varb1 varb3 covb1b3

*     ****************************************************************  *;
*       Calculate data necessary for top marginal effect plot.          *;
*     ****************************************************************  *;

*The first line is important and generates a variable MVZ that takes on 
*all of the values of the modifying variable Z, changenonslav, for which 
*you want to calculate the marginal effect of X; the more values, the

*to start at 0 and go up 0.01= *gen MVZ=((_n-1)/100)
gen MVZ=((_n-.0001)/10000)
*drop MVZ
*The second line determines the last value of the modifying variable Z for which 
*you want to calculate the marginal effect of X. Given how it has been constructed, 
*MVZ = 13 when _n = 1501. Since the highest observed value of Z in my data is 13, 
*I replace MVZ as ., or missing, for all values of _n>1501.
replace  MVZ=. if _n>10000

*The third line calculates this marginal effect of X for all values MVZ, 
*so long as MVZ is less than 13.
gen conbx=b1+b3*MVZ if _n<10001

*The fourth line calculates the standard error for the marginal effect of X 
*for all values MVZ, so long as MVZ is less than 13
gen consx=sqrt(varb1+varb3*(MVZ^2)+2*covb1b3*MVZ) if _n<10001

*The last three lines use the standard errors and marginal effects to create
*two-tailed 95% confidence intervals.
*These lines create the necessary information if you want to include a rug plot 
*in the marginal effect plot

gen ax=1.96*consx

gen upperx=conbx+ax

gen lowerx=conbx-ax

*     ****************************************************************  *;
*       Construct stuff to produce the rug plot.  Need to create an     *;
*       offset position for the pipe marker, which will depend on       *;
*       where the histogram sits on the y-axis. This will requie some   *;
*       trial and error.                                                *;
*     ****************************************************************  *;

*rug plot (have to play around with first line)
*gen where=-1
*drop where

*gen pipe = "|"

*egen tag_nonslav = tag(elecdm)

*     ****************************************************************  *;
*       Construct variable to produce y=0 line.                         *;
*     ****************************************************************  *;

*This will create a horizontal line at 0 on the vertical axis.
gen yline=0

*     ****************************************************************  *;
*     ****************************************************************  *;
*       Produce marginal effect plot for X.                             *;
*     ****************************************************************  *;
*     ****************************************************************  *;

graph twoway hist polyarchy if used==1, width(0.01) percent color(gs14) yaxis(2) ///
	    ||   line conbx   MVZ, clpattern(solid) clwidth(medium) clcolor(black) yaxis(1) ///
        ||   line upperx  MVZ, clpattern(dash) clwidth(thin) clcolor(black) ///
        ||   line lowerx  MVZ, clpattern(dash) clwidth(thin) clcolor(black) ///
        ||   line yline  MVZ,  clwidth(thin) clcolor(black) clpattern(solid) ///
	    ||   , ///
             xlabel(0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1, nogrid labsize(2)) ///
		     ylabel(-4 -3 -2 -1  0 1 2 , axis(1) nogrid labsize(2)) ///
		     ylabel(0 1 2 3 4  , axis(2) nogrid labsize(2)) ///
	         yscale(noline alt) ///
		     yscale(noline alt axis(2)) ///	
             xscale(noline) ///
             legend(off) ///
			 subtitle(" ") ///
             xtitle("Electoral Democracy (V-Dem)" , size(2.5)  ) ///
			 ytitle("Marginal Effect of Judicial Independence" , axis(1) size(2.5)) ///
             ytitle("Percentage of Observations" , axis(2) size(2.5) orientation(rvertical)) ///
             xsca(titlegap(2)) ///
             ysca(titlegap(2)) ///
			 scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white))
